/* "Efficiency and water use: Dynamic effects of irrigation technology adoption"
by Micah Cameron-Harp and Nathan Hendricks

Code written by Micah Cameron-Harp
May 16th, 2024

This do file reads the TWFE, Callaway & Sant'Anna ATT, and 
	de Chaisemartin & D'Haultfoueille (2021) dynamic treatment effect
	estimates and generates figures D.4 and D.5.
*/

/* NOTE - Global macros for directories are created in the appendices_results.do 
	file which calls this individual .do file */

********************************************************************************
/***************************** Flood to CP/LEPA *******************************/
********************************************************************************
*Load the twfe estimates, store the mean of the dependent variables
import excel "${dr_temp}\twfe.xlsx", firstrow sheet(floodtocplepa_91_15_nyt) clear	
		sum panel_mean_dep_var if dep_var=="af_used"
		local mean_af_used = r(mean)
		sum panel_mean_dep_var if dep_var=="acres_irr"
		local mean_acres_irr = r(mean)
		sum panel_mean_dep_var if dep_var=="depth_applied"
		local mean_depth_applied = r(mean)
/* Import the flood to cp or lepa event study coefficients */
use "${dr_temp}\didl_flood_cplepa_af_used.dta", clear
	gen dep_var = "af_used"
	gen panel_mean_dep_var=`mean_af_used'
	*Append the other two dependent variables 
	append using "${dr_temp}\didl_flood_cplepa_acres_irr.dta", force
	replace dep_var = "acres_irr" if missing(dep_var)
	replace panel_mean_dep_var=`mean_acres_irr' if missing(panel_mean_dep_var)
	append using "${dr_temp}\didl_flood_cplepa_depth_applied.dta", force
	replace dep_var = "depth_applied" if missing(dep_var)
	replace panel_mean_dep_var=`mean_depth_applied' if missing(panel_mean_dep_var)
	drop if time_to_treatment==.
	gen estimator = "C&D"
	rename (time_to_treatment treatment_effect se_treatment_effect treatment_effect_upper_95CI treatment_effect_lower_95CI) ///
		(effect_estimated estimate se_estimate ub_estimate lb_estimate)
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep effect_estimated dep_var estimator scaled_estimate scaled_lb scaled_ub estimate lb_estimate ub_estimate

	*Append the TWFE and CS results 
	append using "${dr_temp}\floodtocplepa_es_twfe.dta"
	append using "${dr_temp}\floodtocplepa_es_cs.dta"
 
*Sort 
sort dep_var effect_estimated

*Now offset the location of scaled_estimates for C&S so they don't overlap
replace effect_estimated = effect_estimated + 0.5 if estimator=="C&D"
replace effect_estimated = effect_estimated + 0.25 if estimator=="c&s"

*Replace extreme ub and lb estimates with -50 and +50 to make graphs readable 
gen extreme = scaled_ub > 40 | scaled_lb < -40
replace scaled_ub = 40 if scaled_ub > 40
replace scaled_lb = -40 if scaled_lb < -40


/* Graph C&S, C&D, and TWFE together */
*AF_USED_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="af_used" & extreme!=1, msize(.2) lwidth(.1) lcolor(orange%60)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="af_used", lwidth(.1) lcolor(orange%60) msize(.4) mcolor(orange%60)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="af_used", msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%60)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & dep_var=="af_used" & extreme!=1, msize(.2) lwidth(.1) lcolor(purple%80)) /// 
		(connected scaled_estimate effect_estimated if estimator=="C&D" & dep_var=="af_used", lwidth(.1) lcolor(purple%80) msize(.4) mcolor(purple%80) msymbol(D)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & extreme==1 & dep_var=="af_used", msize(0) lpattern("shortdash") lwidth(.1) lcolor(purple%80)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="af_used" & extreme!=1,  msize(.2) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="af_used", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="af_used", msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Flood to center pivot", size(small)) ///
			subtitle("Acre-feet withdrawn", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-24 -20(5)20 23, nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(8 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(5 "{it:{&delta}{superscript:CD}(l)} from de Chaisemartin & D'Haultfoueille (2022)")  ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
			graphregion(color(white) margin(zero)) name(es_af_used, replace)

*ACRES_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="acres_irr" & extreme!=1, msize(.2) lwidth(.1) lcolor(orange%60)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="acres_irr", lwidth(.1) lcolor(orange%60) msize(.4) mcolor(orange%60)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="acres_irr", msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%60)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & dep_var=="acres_irr" & extreme!=1, msize(.2) lwidth(.1) lcolor(purple%80)) /// 
		(connected scaled_estimate effect_estimated if estimator=="C&D" & dep_var=="acres_irr", lwidth(.1) lcolor(purple%80) msize(.4) mcolor(purple%80) msymbol(D)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & extreme==1 & dep_var=="acres_irr", msize(0) lpattern("shortdash") lwidth(.1) lcolor(purple%80)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="acres_irr" & extreme!=1,  msize(.2) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="acres_irr", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="acres_irr", msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Flood to center pivot", size(small)) ///
			subtitle("Acres irrigated", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-24 -20(5)20 23, nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(8 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(5 "{it:{&delta}{superscript:CD}(l)} from de Chaisemartin & D'Haultfoueille (2022)")  ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_acres_irr, replace)

*DEPTH APPLIED
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="depth_applied" & extreme!=1, msize(.2) lwidth(.1) lcolor(orange%60)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="depth_applied", lwidth(.1) lcolor(orange%60) msize(.4) mcolor(orange%60)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="depth_applied", msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%60)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & dep_var=="depth_applied" & extreme!=1, msize(.2) lwidth(.1) lcolor(purple%80)) /// 
		(connected scaled_estimate effect_estimated if estimator=="C&D" & dep_var=="depth_applied", lwidth(.1) lcolor(purple%80) msize(.4) mcolor(purple%80) msymbol(D)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & extreme==1 & dep_var=="depth_applied", msize(0) lpattern("shortdash") lwidth(.1) lcolor(purple%80)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="depth_applied" & extreme!=1,  msize(.2) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="depth_applied", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="depth_applied", msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Flood to center pivot", size(small)) ///
			subtitle("Depth applied", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-24 -20(5)20 23, nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(8 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(5 "{it:{&delta}{superscript:CD}(l)} from de Chaisemartin & D'Haultfoueille (2022)")  ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_depth_applied, replace)
			
/* Combine the graphs */
grc1leg2 es_af_used es_acres_irr es_depth_applied, ///
	rows(3) cols(1) legendfrom(es_af_used) ///
	maintotoptitle maintitlefrom(es_af_used) ///
	xtob1title xtitlefrom(es_af_used) ///
	ytol1title ytitlefrom(es_af_used) ///
	labsize(vsmall) graphregion(color(white)) xsize(6.5) ysize(8) iscale(.8) ring(1)
graph export "${dr_output_app}/figureD4.tif", width(6500) height(7900) replace

********************************************************************************
/***************************** CP to LEPA *******************************/	
********************************************************************************
*Load the twfe estimates, store the mean of the dependent variables
import excel "${dr_temp}\twfe.xlsx", firstrow sheet(cptolepa_91_19) clear
		sum panel_mean_dep_var if dep_var=="af_used"
		local mean_af_used = r(mean)
		sum panel_mean_dep_var if dep_var=="acres_irr"
		local mean_acres_irr = r(mean)
		sum panel_mean_dep_var if dep_var=="depth_applied"
		local mean_depth_applied = r(mean)
/* Import the cp to lepa event study coefficients */
	use "${dr_temp}\didl_cp_lepa_af_used.dta", clear
	gen dep_var = "af_used"
	gen panel_mean_dep_var=`mean_af_used'
	*Append the other two dependent variables 
	append using "${dr_temp}\didl_cp_lepa_acres_irr.dta", force
	replace dep_var = "acres_irr" if missing(dep_var)
	replace panel_mean_dep_var=`mean_acres_irr' if missing(panel_mean_dep_var)
	append using "${dr_temp}\didl_cp_lepa_depth_applied.dta", force
	replace dep_var = "depth_applied" if missing(dep_var)
	replace panel_mean_dep_var=`mean_depth_applied' if missing(panel_mean_dep_var)
	drop if time_to_treatment==.
	gen estimator = "C&D"
	rename (time_to_treatment treatment_effect se_treatment_effect treatment_effect_upper_95CI treatment_effect_lower_95CI) ///
		(effect_estimated estimate se_estimate ub_estimate lb_estimate)
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep effect_estimated dep_var estimator scaled_estimate scaled_lb scaled_ub estimate lb_estimate ub_estimate

	*Append the TWFE and CS results 
	append using "${dr_temp}\cptolepa_es_twfe.dta"
	append using "${dr_temp}\cptolepa_es_cs.dta"
 
*Sort 
sort dep_var effect_estimated

*Now offset the location of scaled_estimates for C&S so they don't overlap
replace effect_estimated = effect_estimated + 0.5 if estimator=="C&D"
replace effect_estimated = effect_estimated + 0.25 if estimator=="c&s"

*Replace extreme ub and lb estimates with -50 and +50 to make graphs readable 
gen extreme = scaled_ub > 40 | scaled_lb < -40
replace scaled_ub = 40 if scaled_ub > 40
replace scaled_lb = -40 if scaled_lb < -40

/* Graph C&S together with C&D */
*AF_USED_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="af_used" & extreme!=1, msize(.2) lwidth(.1) lcolor(orange%60)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="af_used", lwidth(.1) lcolor(orange%60) msize(.4) mcolor(orange%60)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="af_used", msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%60)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & dep_var=="af_used" & extreme!=1, msize(.2) lwidth(.1) lcolor(purple%80)) /// 
		(connected scaled_estimate effect_estimated if estimator=="C&D" & dep_var=="af_used", lwidth(.1) lcolor(purple%80) msize(.4) mcolor(purple%80) msymbol(D)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & extreme==1 & dep_var=="af_used", msize(0) lpattern("shortdash") lwidth(.1) lcolor(purple%80)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="af_used" & extreme!=1,  msize(.2) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="af_used", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="af_used", msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Traditional center pivot to LEPA", size(small)) ///
			subtitle("Acre-feet withdrawn", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-28 -20(5)20 27,  nogrid) ///
			xscale(range(-28.5 27.5)) ///
			legend(order(8 5 2) label(8 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(5 "{it:{&delta}{superscript:CD}(l)} from de Chaisemartin & D'Haultfoueille (2022)")  ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
			graphregion(color(white) margin(zero)) name(es_af_used, replace)

*ACRES_IRR
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="acres_irr" & extreme!=1, msize(.2) lwidth(.1) lcolor(orange%60)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="acres_irr", lwidth(.1) lcolor(orange%60) msize(.4) mcolor(orange%60)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="acres_irr", msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%60)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & dep_var=="acres_irr" & extreme!=1, msize(.2) lwidth(.1) lcolor(purple%80)) /// 
		(connected scaled_estimate effect_estimated if estimator=="C&D" & dep_var=="acres_irr", lwidth(.1) lcolor(purple%80) msize(.4) mcolor(purple%80) msymbol(D)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & extreme==1 & dep_var=="acres_irr", msize(0) lpattern("shortdash") lwidth(.1) lcolor(purple%80)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="acres_irr" & extreme!=1,  msize(.2) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="acres_irr", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="acres_irr", msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Traditional center pivot to LEPA", size(small)) ///
			subtitle("Acres irrigated", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-28 -20(5)20 27,  nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(8 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(5 "{it:{&delta}{superscript:CD}(l)} from de Chaisemartin & D'Haultfoueille (2022)")  ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_acres_irr, replace)

*DEPTH APPLIED
	tw 	(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & dep_var=="depth_applied" & extreme!=1, msize(.2) lwidth(.1) lcolor(orange%60)) /// 
		(connected scaled_estimate effect_estimated if estimator=="twfe" & dep_var=="depth_applied", lwidth(.1) lcolor(orange%60) msize(.4) mcolor(orange%60)  msymbol(T)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="twfe" & extreme==1 & dep_var=="depth_applied", msize(0) lpattern("shortdash") lwidth(.1) lcolor(orange%60)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & dep_var=="depth_applied" & extreme!=1, msize(.2) lwidth(.1) lcolor(purple%80)) /// 
		(connected scaled_estimate effect_estimated if estimator=="C&D" & dep_var=="depth_applied", lwidth(.1) lcolor(purple%80) msize(.4) mcolor(purple%80) msymbol(D)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="C&D" & extreme==1 & dep_var=="depth_applied", msize(0) lpattern("shortdash") lwidth(.1) lcolor(purple%80)) /// 
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & dep_var=="depth_applied" & extreme!=1,  msize(.2) lwidth(.1) lcolor(blue%100)) /// 
		(connected scaled_estimate effect_estimated if estimator=="c&s" & dep_var=="depth_applied", lwidth(.1) lcolor(blue%100) msize(.4) mcolor(blue%100) msymbol(O)) ///
		(rcap scaled_ub scaled_lb effect_estimated if estimator=="c&s" & extreme==1 & dep_var=="depth_applied", msize(0) lpattern("shortdash") lwidth(.1) lcolor(blue%100)), /// 
			title("        Traditional center pivot to LEPA", size(small)) ///
			subtitle("Depth applied", size(small)) ///
			xtitle("Time relative to adoption (years)", size(vsmall)) ///
			yline(0, lpattern(solid)) ///
			ytitle("Percent change", size(small)) ///
			ylab(-40(20)40, ang(45) nogrid) ///
			xlab(-28 -20(5)20 27,  nogrid) ///
			xscale(range(-24 24)) ///
			legend(order(8 5 2) label(8 "{it:{&delta}{superscript:CS}(l)} from Callaway & Sant'Anna (2020)") ///
							  label(5 "{it:{&delta}{superscript:CD}(l)} from de Chaisemartin & D'Haultfoueille (2022)")  ///
							  label(2 "{it:{&beta}{superscript:TWFE}(l)} TWFE event study")  ///
								cols(1) size(vsmall)) ///
								graphregion(color(white) margin(zero)) name(es_depth_applied, replace)
			
/* Combine the graphs */
grc1leg2 es_af_used es_acres_irr es_depth_applied, ///
	rows(3) cols(1) legendfrom(es_af_used) ///
	maintotoptitle maintitlefrom(es_af_used) ///
	xtob1title xtitlefrom(es_af_used) ///
	ytol1title ytitlefrom(es_af_used) ///
	labsize(vsmall) graphregion(color(white)) xsize(6.5) ysize(8) iscale(.8) ring(1)
graph export "${dr_output_app}/figureD5.tif", width(6500) height(7900) replace